function [out,factorSim,yieldSim] = simulator(setupStep3,theta11Step3,theta2Step3,matSelectSixTh,numSim)
nx = setupStep3.nx;
[h0,hx,sigmaVar]   = unfoldTheta2(theta2Step3,nx);

%% Simulation
factorSim      = zeros(nx,numSim);
factorSim(:,1) = (eye(nx)-hx)\h0; %inv(eye(nx)-hx)*h0
randn('seed',1)
epsilon        = randn(nx,numSim);
epsilon        = (epsilon-repmat(mean(epsilon,2),1,numSim))./repmat(std(epsilon,0,2),1,numSim);
for s=2:numSim
    factorSim(:,s) = h0 + hx*factorSim(:,s-1) + sigmaVar*epsilon(:,s);
end

% Yields and premia for all yields
[yHatTrans,termPremSim,fHatSim,fwdPremSim] = getTermPremium(theta11Step3,theta2Step3,setupStep3,factorSim,1);
yieldSim = yHatTrans';

% Empirical moments
meanYield = mean(yHatTrans,1);
stdYield = std(yHatTrans,0,1);
LPYI = CampbellSchillerRegression(yHatTrans,zeros(size(yHatTrans)),zeros(size(yHatTrans)),1);
LPYI.CSBetta       = LPYI.CSBetta(:,matSelectSixTh);
LPYI.CSBettaAdj    = LPYI.CSBettaAdj(:,matSelectSixTh);
LPYI.CSBetta_se    = LPYI.CSBetta_se(:,matSelectSixTh);
LPYI.CSBettaAdj_se = LPYI.CSBettaAdj_se(:,matSelectSixTh);
%LPYII = CampbellSchillerRegression(yHatTrans,termPremSim,fwdPremSim,1);
[skew,kurt]        = skewnessAndKurtosis(yHatTrans(2:end,:)-yHatTrans(1:end-1,:));
yieldReduced       = yieldSim(matSelectSixTh,:);
slope.mean         = mean(yieldReduced-repmat(yieldReduced(1,:),size(yieldReduced,1),1),2)';
slope.std          = std(yieldReduced-repmat(yieldReduced(1,:),size(yieldReduced,1),1),0,2)';
%% Saving output
out.meanYield   = meanYield(matSelectSixTh);
out.stdYield    = stdYield(matSelectSixTh);
out.skewYield   = skew(matSelectSixTh);
out.kurtYield   = kurt(matSelectSixTh);
out.slope       = slope;
out.LPYI        = LPYI;
%out.LPYII       = LPYII;
out.meanTP      = mean(termPremSim(:,matSelectSixTh),1);
out.stdTP       = std(termPremSim(:,matSelectSixTh),0,1);
out.numSim      = numSim;
out.matSelect   = matSelectSixTh;
%out.factorSim   = factorSim;
%out.yieldSim    = yieldSim;

end